mrgsolveR package for simulation from ODE-based models
C++ inside model specification
formatODEPACK / DLSODA
(FORTRAN)
RID, \(\eta\), \(\varepsilon\)F, ALAG, SS
etc, handled under the hoodR, C++,
Rcpp, boost, RcppArmadilloR is it’s natural habitatmrgsolve started as QSP modeling toolR front end to deSolveC++ interface to DLSODARcpp - vectors, matrices, functions, environments,
random numbersboost - numerical tools in C++C++ code (functions, data structures,
classes)SBML to mrgsolve using
R bindings to libSBMLGitHub site: https://github.com/metrumresearchgroup/mrgsolve
Discussion: https://github.com/metrumresearchgroup/mrgsolve/discussions
Issues: https://github.com/metrumresearchgroup/mrgsolve/issues
mrgsolve website: https://mrgsolve.org
User Guide: https://mrgsolve.org/user_guide
Vignettes: https://mrgsolve.org/vignettes
mod %>% ev(amt = 100, ii = 24, addl = 3) %>% mrgsim() %>% plot()
mod %>% ev(amt = 100, ii = 24, addl = 3) %>% mrgsim() %>% plot()
that come from?event in this exampleplot, mutate, as_tibble etc
…
model %>%intervention %>%Go! %>%take-a-look
model %>% intervention %>%options %>% Go! %>% ...
model %>% intervention %>%population %>% Go! %>% ...
model %>%data-set %>% Go! %>% ...
- where
data-set =intervention + population
- For now, let’s get this part down
model %>%intervention %>%Go! %>%take-a-look
%>% ?What happens first in this operation?
mean(sqrt(seq(4)))
mean(sqrt(seq(4)))
4 %>% seq() %>% sqrt() %>% mean()
Better.
4 %>%
seq(.) %>%
sqrt(.) %>%
mean(., na.rm = TRUE)
mod %>% some_intervention() %>% simulate() %>% post_process()
%>% ...
- I (almost) always call the model object
mod in the documention / examples
- All the information about the model we need to know to run the simulation
Distinct from the intervention, the population, the summary
- But the model
does know about output time, random effects
mod
----------------- source: pk1.cpp -----------------
project: /Users/ahmede/Li...gsolve/models
shared object: pk1-so-853b215d7005
time: start: 0 end: 192 delta: 0.2
add: <none>
compartments: EV CENT [2]
parameters: CL V KA [3]
captures: CP [1]
omega: 0x0
sigma: 0x0
solver: atol: 1e-08 rtol: 1e-08 maxsteps: 20k
------------------------------------------------------
param(mod)
Model parameters (N=3):
name value . name value
CL 1 | V 20
KA 1 | . .
- Parameters get a name
- Names and number of parameters gets fixed at compile time
- All parameters have a value
- Value can be modified after compile time
init(mod)
Model initial conditions (N=2):
name value . name value
CENT (2) 0 | EV (1) 0
- Every compartment gets a name
- Every compartment gets an initial condition
mod <- mread("simple", "model")
mod <- mread("<model-name>", "<project-directory>")
First time to read
mod <- mread_cache("simple", here("day1/model"))
. Building simple ... done.
Next time to read
mod <- mread_cache("simple", here("day1/model"))
. Loading model from cache.
soloc)It can be helpful to set
options(mrgsolve.soloc = "my_build_directory")
Equivalent to
mod <- mread("simple", here("day1/model"), soloc = "my_build_directory")
mod <- mread("<first-argument>", "<second-argument>")
mod <- mread("<first-argument>", modlib())
modlib()
. [1] "/Users/ahmede/Library/Caches/org.R-project.R/R/renv/cache/v5/macos/R-4.4/aarch64-apple-darwin23.4.0/mrgsolve/1.5.1/689a814a168f375aa0ebfeeb3ee78394/mrgsolve/models"
mod <- mread("effect", modlib())
mod
.
.
. ---------------- source: effect.cpp ----------------
.
. project: /Users/ahmede/Li...gsolve/models
. shared object: effect-so-b1e633a5a477
.
. time: start: 0 end: 36 delta: 0.1
. add: <none>
.
. compartments: GUT CENT PERIPH CE [4]
. parameters: VC KA K10 K12 K21 E0 EMAX EC50 KEO [9]
. captures: EFFECT CP [2]
. omega: 0x0
. sigma: 0x0
.
. solver: atol: 1e-08 rtol: 1e-08 maxsteps: 20k
. ------------------------------------------------------
Open R script in Rstudio then …
library(mrgsolve)
mod <- mread("pbpk", modlib())
mod <- modlib("pbpk")
File name:
Section name:
model %>% %>% Go! %>% take-a-look
Event object = quick / easy way to implement dose or other intervention
e <- ev(amt = 100)
e
. Events:
. time amt cmt evid
. 1 0 100 1 1
- Defaults:
time,evid,cmt
mod %>% ev(amt = 100) %>% mrgsim()
e <- ev(amt = 100)
mod %>% ev(e) %>% mrgsim()
mod %>% mrgsim(events = e)
ev(...)evidevid 1, with
rate==0)evid 1, with
rate > 0)evid 2)
evid 3)evid 4)evid 8)What’s going to happen?
e1 <- ev(amt = 200)
e2 <- ev(amt = 100, time = 24, ii = 24, addl = 4)
c(e1, e2)
What’s going to happen?
Combine two events
e1 <- ev(amt = 200)
e2 <- ev(amt = 100, time = 24, ii = 24, addl = 4)
c(e1, e2)
. Events:
. time amt cmt evid ii addl
. 1 0 200 1 1 0 0
. 2 24 100 1 1 24 4
What’s going to happen?
e1 <- ev(amt = 200, ii = 12, addl = 2)
e2 <- ev(amt = 100, ii = 24, addl = 4)
seq(e1, e2)
What’s going to happen?
Put two events in a sequence
e1 <- ev(amt = 200, ii = 12, addl = 2)
e2 <- ev(amt = 100, ii = 24, addl = 4)
seq(e1, e2)
. Events:
. time amt ii addl cmt evid
. 1 0 200 12 2 1 1
. 2 36 100 24 4 1 1
What is going to happen?
e1 <- ev(amt = 200)
e2 <- ev(amt = 100, ii = 24, addl = 4)
seq(e1, wait = 36, e2)
What is going to happen?
Wait before starting the next part of the regimen
e1 <- ev(amt = 200)
e2 <- ev(amt = 100, ii = 24, addl = 4)
seq(e1, wait = 36, e2)
. Events:
. time amt ii addl cmt evid
. 1 0 200 0 0 1 1
. 2 36 100 24 4 1 1
day1/examples/pbpk-rifampin-ddi.Ras.data.frame(e1)
. time amt cmt evid
. 1 0 200 1 1
We will use a
Event objects are convenient
model %>% intervention %>%
e <- ev(amt = 100)
mod %>% ev(e) %>% mrgsim()
. Model: effect
. Dim: 362 x 8
. Time: 0 to 36
. ID: 1
. ID time GUT CENT PERIPH CE EFFECT CP
. 1: 1 0.0 0.00 0.000 0.0000 0.0000 157.0 0.000
. 2: 1 0.0 100.00 0.000 0.0000 0.0000 157.0 0.000
. 3: 1 0.1 91.21 8.443 0.1552 0.2225 155.7 3.460
. 4: 1 0.2 83.19 15.502 0.5817 0.8048 152.8 6.353
. 5: 1 0.3 75.88 21.354 1.2272 1.6382 149.6 8.752
. 6: 1 0.4 69.21 26.156 2.0463 2.6356 146.6 10.719
. 7: 1 0.5 63.13 30.046 3.0001 3.7278 144.1 12.314
. 8: 1 0.6 57.58 33.146 4.0553 4.8607 142.2 13.584
mrgsim(mod, events = e)
mod %>% ev(amt = 100) %>% mrgsim() %>% plot(CP~time)
mod %>% ev(amt = 100) %>% mrgsim(end = 48) %>% plot(CP~time)
mod %>% ev(amt = 100) %>% mrgsim(end = 48, delta = 0.1) %>% plot(CP~time)
0stime(mod)
. [1] 0 6 12 18
%>% %>% %>%
mread or
mread_cachemread("<model-name>", modlib())param(mod)init(mod)see(mod)ev(...): amt, cmt,
time, ii, addl,
ratemodel %>% intervention %>% Go! %>%
mod <- mread_cache("effect", modlib())
mod %>% ev(amt = 100) %>% mrgsim() %>% plot()
mod %>% ev(amt = 100) %>% mrgsim() %>%
plot(CP + EFFECT~., col = "firebrick")
. Warning: In subset.data.frame(as.data.frame(x), ...) :
. extra argument 'col' will be disregarded
mod %>% mrgsim() %>% mutate(arm = 1)
. # A tibble: 361 × 9
. ID time GUT CENT PERIPH CE EFFECT CP arm
. <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
. 1 1 0 0 0 0 0 157 0 1
. 2 1 0.1 0 0 0 0 157 0 1
. 3 1 0.2 0 0 0 0 157 0 1
. 4 1 0.3 0 0 0 0 157 0 1
. 5 1 0.4 0 0 0 0 157 0 1
. 6 1 0.5 0 0 0 0 157 0 1
. 7 1 0.6 0 0 0 0 157 0 1
. 8 1 0.7 0 0 0 0 157 0 1
. 9 1 0.8 0 0 0 0 157 0 1
. 10 1 0.9 0 0 0 0 157 0 1
. # ℹ 351 more rows
mod %>% mrgsim() %>% head(n = 2)
. ID time GUT CENT PERIPH CE EFFECT CP
. 1 1 0.0 0 0 0 0 157 0
. 2 1 0.1 0 0 0 0 157 0
mod %>% Req(CP) %>% mrgsim() %>% head(n = 2)
. ID time CP
. 1 1 0.0 0
. 2 1 0.1 0
mod %>% ev(amt=1) %>% mrgsim() %>% head(n = 2)
. ID time GUT CENT PERIPH CE EFFECT CP
. 1 1 0 0 0 0 0 157 0
. 2 1 0 1 0 0 0 157 0
mod %>% ev(amt=1) %>% obsonly() %>% mrgsim() %>% head(n = 2)
. ID time GUT CENT PERIPH CE EFFECT CP
. 1 1 0.0 0.0000000 0.00000000 0.000000000 0.000000000 157.0000 0.00000000
. 2 1 0.1 0.9121052 0.08443153 0.001551552 0.002224565 156.9866 0.03460309
On the fly
mod %>% update(end = 72) %>% mrgsim()
Persistent
mod <- update(mod, end = 72)
mrgsim will call update for you (on the fly)
mod %>% mrgsim(end = 72)
start, end, delta,
addatol, rtolhmax, maxsteps, mxhnil,
ixpr$OMEGA, $SIGMAtscale (rescale the output time)digitsmod %>% update(rtol = 1E-12) %>% ...
mod %>% mrgsim(rtol = 1E-12) %>% ...
mod <- param(mod, CL = 1.5)
%>% %>% %>%
model %>% intervention %>%population %>% Go! %>% take-a-look
mod %>% ev(amt = 100, ii = 24, addl = 3) %>%idata_set( pop ) %>% mrgsim() %>% plot()
idata_set takes in individual-level datahead(pop, n = 3)
. ID CL VC KA KOUT IC50 FOO
. 1 1 1.05 47.8 0.8390 2.45 1.28 4
. 2 2 0.73 30.1 0.0684 2.51 1.84 6
. 3 3 2.82 23.8 0.1180 3.88 2.48 5
length(unique(pop$ID))
. [1] 10
head(pop, n = 3)
. ID CL VC KA KOUT IC50 FOO
. 1 1 1.05 47.8 0.8390 2.45 1.28 4
. 2 2 0.73 30.1 0.0684 2.51 1.84 6
. 3 3 2.82 23.8 0.1180 3.88 2.48 5
param(mod)
.
. Model parameters (N=3):
. name value . name value
. CL 1 | V 20
. KA 1 | . .
idata?Batches of simulations or sensitivity analyses
idata <- expand.idata(CL = seq(0.5, 1.5, 0.25))
idata
. ID CL
. 1 1 0.50
. 2 2 0.75
. 3 3 1.00
. 4 4 1.25
. 5 5 1.50
+
mod %>%
idata_set(idata) %>%
ev(amt = 100, ii = 24, addl = 2) %>%
mrgsim(end = 120) %>% plot(log(CP) ~ .)
Description: sensitivity analysis in a PBPK model
File name:
data_set is the dosing equivalent to
idata_setdata <- expand.ev(amt = c(100, 300, 1000), ii = 24, addl = 3)
head(data)
. ID time amt ii addl cmt evid
. 1 1 0 100 24 3 1 1
. 2 2 0 300 24 3 1 1
. 3 3 0 1000 24 3 1 1
data_set configurationmod %>%
data_set(data) %>%
mrgsim(end = 120) %>% plot(log(CP) ~ .)
data <- expand.ev(
amt = c(100,300,1000),
ii = 24, addl = 3,
CL = seq(0.5,1.5, 0.5)
)
head(data)
. ID time amt ii addl cmt evid CL
. 1 1 0 100 24 3 1 1 0.5
. 2 2 0 300 24 3 1 1 0.5
. 3 3 0 1000 24 3 1 1 0.5
. 4 4 0 100 24 3 1 1 1.0
. 5 5 0 300 24 3 1 1 1.0
. 6 6 0 1000 24 3 1 1 1.0
data <- mutate(data, dose = amt)
mod %>%
data_set(data) %>%
mrgsim(carry.out = "dose", end = 120) %>%
plot(log(CP)~time|factor(dose), group = ID, scales = "same")
Copy from the input data to the output data
mod %>% carry_out(dose) %>% data_set(data, dose==300) %>% mrgsim()
. Model: pk1
. Dim: 2886 x 4
. Time: 0 to 192
. ID: 3
. ID time dose CP
. 1: 2 0.0 300 0.000
. 2: 2 0.0 300 0.000
. 3: 2 0.2 300 2.712
. 4: 2 0.4 300 4.919
. 5: 2 0.6 300 6.712
. 6: 2 0.8 300 8.167
. 7: 2 1.0 300 9.345
. 8: 2 1.2 300 10.296
head(data, n = 3)
. ID time amt ii addl cmt evid CL dose
. 1 1 0 100 24 3 1 1 0.5 100
. 2 2 0 300 24 3 1 1 0.5 300
. 3 3 0 1000 24 3 1 1 0.5 1000
ev(amt = 100, ii = 24, addl = 3)
. Events:
. time amt ii addl cmt evid
. 1 0 100 24 3 1 1
mod + ev = ?mod + idata_set + ev = ?mod + data_set = ?
head(data)
. ID time amt ii addl cmt evid CL dose
. 1 1 0 100 24 3 1 1 0.5 100
. 2 2 0 300 24 3 1 1 0.5 300
. 3 3 0 1000 24 3 1 1 0.5 1000
. 4 4 0 100 24 3 1 1 1.0 100
. 5 5 0 300 24 3 1 1 1.0 300
. 6 6 0 1000 24 3 1 1 1.0 1000
Recall that data sets are just plain old data.frames in R. Feel free to code these up in any way that is convenient for you.
In our experience, the following helper functions cover many (not every) common needs for building the data sets.
expand.evev_repas_data_setev_assignexpand.evLike expand.grid
expand.ev(ID = 1:2, amt = c(100,200))
. ID time amt cmt evid
. 1 1 0 100 1 1
. 2 2 0 100 1 1
. 3 3 0 200 1 1
. 4 4 0 200 1 1
ev_repe <- ev(amt = 100, ii = 12, addl = 14)
ev_rep(e, ID = seq(5))
. ID time amt ii addl cmt evid
. 1 1 0 100 12 14 1 1
. 2 2 0 100 12 14 1 1
. 3 3 0 100 12 14 1 1
. 4 4 0 100 12 14 1 1
. 5 5 0 100 12 14 1 1
e <- ev(amt = 100, ii = 12, addl = 14, ID = seq(5))
ev_repe <- seq(ev(amt = 100), wait = 36, ev(amt = 50, ii = 24, addl = 4))
e
. Events:
. time amt ii addl cmt evid
. 1 0 100 0 0 1 1
. 2 36 50 24 4 1 1
ev_rep(e, ID = seq(2))
. ID time amt ii addl cmt evid
. 1 1 0 100 0 0 1 1
. 2 1 36 50 24 4 1 1
. 3 2 0 100 0 0 1 1
. 4 2 36 50 24 4 1 1
as_data_setdata <- as_data_set(
ev(amt = 100, ii = 12, addl = 19, ID = 1:2),
ev(amt = 200, ii = 24, addl = 9, ID = 1:3),
ev(amt = 150, ii = 24, addl = 9, ID = 1:4)
)
. ID time amt ii addl cmt evid
. 1 1 0 100 12 19 1 1
. 2 2 0 100 12 19 1 1
. 3 3 0 200 24 9 1 1
. 4 4 0 200 24 9 1 1
. 5 5 0 200 24 9 1 1
. 6 6 0 150 24 9 1 1
. 7 7 0 150 24 9 1 1
. 8 8 0 150 24 9 1 1
. 9 9 0 150 24 9 1 1
ev_assignpop <- expand.idata(GROUP = c(1,2), ID = 1:3)
head(pop)
. ID GROUP
. 1 1 1
. 2 2 2
. 3 3 1
. 4 4 2
. 5 5 1
. 6 6 2
e1 <- ev(amt = 100, ii = 24, addl = 9)
e2 <- mutate(e1, amt = 200)
data <- ev_assign(
l = list(e1,e2),
idata = pop,
evgroup = "GROUP",
join = TRUE
)
head(data)
. time amt ii addl cmt evid ID GROUP
. 1 0 100 24 9 1 1 1 1
. 2 0 200 24 9 1 1 2 2
. 3 0 100 24 9 1 1 3 1
. 4 0 200 24 9 1 1 4 2
. 5 0 100 24 9 1 1 5 1
. 6 0 200 24 9 1 1 6 2
Description: simulate GCSF PK/PD profiles for weight-based dosing
File name:
Description: Population PK of azithromycin
File name:
Description: simulate PK profiles using empirical Bayes estimates from a meropenem model fit in NONMEM
File name:
CL = 1.25ke = CL/V
- What are the
parameters ?- What are the
compartments ?- What are the
differential equations ?
$PARAM [R] declare and
initialize model parameters$PARAM CL = 1, VC = 20, KA=1.2
name = value format, or <newline>name except for words in
mrgsolve:::reserved()R parsername anywhere in the modelR but read-only in the model specification
file
$FIXED if appropriate$PARAM with
the names in your input data sets$CMT [txt] declare model
compartmentsname and number of compartments in
the modelmrgsolve:::reserved()$INITExample $CMT
$CMT GUT CENT
Example $INIT
$INIT GUT = 0, CENT = 10
$ODE [C++] write differential
equationsdxdt_CMT = ratein - rateout;CMT for compartment amounts, parameters, or other
variables$ODE
dxdt_GUT = -KA*GUT;
dxdt_CENT = KA*GUT - (CL/VC)*CENT;
new variables in your modelC++, which is a
language$MAIN, $ODE, $TABLE,
$GLOBALTo get a numeric variable, use
double
double x = 5.4;
Other types you might use
bool cured = false;
int i = 2;
If in doubt, use double; it’s what you want most of the
time.
$MAIN [C++] covariate model
& initialsFor example
$MAIN
double CL = TVCL*pow(WT/70,0.75);
double VC = TVVC*pow(0.86,SEX);
RESP_0 = KIN/KOUT;
In this example
TVCL, KIN, KOUT,
WT, SEX were declared in
$PARAMRESPC++ expressions and functionsif(a == 2) b = 2;
if(b <= 2) {
c=3;
} else {
c=4;
}
double d = pow(base,exponent);
double e = exp(3);
double f = fabs(-4);
double g = sqrt(5);
double h = log(6);
double i = log10(7);
double j = floor(4.2);
double k = ceil(4.2);
Lots of help on the web http://en.cppreference.com/w/cpp/numeric/math/tgamma
#define)$GLOBAL
#define CP (CENT/VC)
#define g 9.8
$ODE
double INH = CP/(EC50+CP);
dxdt_RESP = KIN*(1-INH) - KOUT*RESP;
$CAPTURE CP
POPULATION simulation$OMEGA and
$SIGMA to define covariance matrices for IIV and RUV,
respectively
ETA(n) and
EPS(m) for realized random effects drawn from
$OMEGA and $SIGMA respectively
ETAs insteadmrgsolve recognizes
ID as a subject indicator for simulating a new
ETA
mrgsolve allows you to write
an error model as a function of EPS(m) and any other
calculated value in your model
$OMEGA and $SIGMA
[txt]$OMEGA and $SIGMA are allowed and
each may be namedmrgsolve will
enforce later on$OMEGA 0.1 0.2 0.3
$OMEGA and $SIGMA@ delimited & must be on different line
from the matrix data$OMEGA @block
0.1 0.0947 0.2
@block)$OMEGA @cor
0.1 0.67 0.2
$OMEGA and $SIGMA$OMEGA and $SIGMA are allowed;
each may be named$OMEGA @name PK
0 0 0
$OMEGA @name PD
0 0
Users are encouraged to add labels
$OMEGA @name PK @labels ECL EVC EKA
0 0 0
@ macros$OMEGA and
$SIGMA@ should appear at the start of the line and the entire
line is reserved for options@cor @name PKTwo forms:
@block means block=TRUE@labels ECL EVC means
labels=c("ECL", "EVC")$PKMODELncmt to 1 or 2depot to TRUE for extravascular dosing
compartmenttrans to change the names of required
parameters$CMT to declare 1 to 3 compartments as
appropriate$PKMODEL cmt = "GUT CENT PERIPH", depot=TRUE
$PARAM CL = 1 , V2 = 30, Q = 8, V3 = 400, KA=1.2
$PKMODEL$PKMODEL ncmt=1, depot=TRUE
$CMT GUT CENT
$PARAM TVCL = 1 , TVV = 30, KA=1.2
$OMEGA @labels ECL EVC
1 1
$MAIN
double CL = TVCL*exp(ECL);
double V = TVV *exp(EVC);
To change the bioavability of doses administered to a
compartment, set F_CMT in $MAIN
To add a lagtime for doses administered to a compartment, set
ALAG_CMT in $MAIN
To set the duration of infusion, set D_CMT in
$MAIN
rate = -2 in your data set or event
objectTo set the infusion rate, set R_CMT in
$MAIN
rate = -1 in your data set or event
objectx, and yzadd <- function(x,y) {
z <- x + y
return(z)
}
add(x = 2, y = 3)
. [1] 5
return() is not strictly necessaryadd <- function(x, y) {
x + y
}
add(2,3)
. [1] 5
add <- function(x, y=3) {
x + y
}
add(2)
. [1] 5
We use the lsa() function from the mrgsim.sa package
library(mrgsim.sa)
?lsa
sens <- lsa(model_object, parameter_names, output_variables)
mod <- modlib("irm1", end = 72, delta = 0.1)
mod %>%
ev(amt = 100) %>%
mrgsim(end = 72, delta = 0.1) %>%
plot()
sens <- lsa(
mod = ev(mod, amt = 100),
par = "CL,V2,Q,KA",
var = "CP"
)
head(sens, n=3)
. # A tibble: 3 × 5
. time dv_name dv_value p_name sens
. <dbl> <chr> <dbl> <chr> <dbl>
. 1 0 CP 0 CL 0
. 2 0 CP 0 CL 0
. 3 0.1 CP 0.472 CL -0.00253
plot(sens)
We’ll look at Sobol’s method
sensitivity::sobol2007